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Abstract. Helioseismic measurements indicate that the solar tachocline is very 
thin, its full thickness not exceeding 4% of the solar radius. The mechanism that 
inhibits differential rotation to propagate from the convective zone to deeper into 
the radiative zone is not known, though several propositions have been made. In 
this paper we demonstrate by numerical models and analytic estimates that the 
tachocline can be confined to its observed thickness by a poloidal magnetic field 
Bp of about one kilogauss, penetrating below the convective zone and oscillating 
with a period of 22 years, if the tachocline region is turbulent with a diffusivity 
of r] ~ 10 10 cm 2 /s (for a turbulent magnetic Prandtl number of unity). We also 
show that a similar confinement may be produced for other pairs of the parameter 
values (Bp, rf). The assumption of the dynamo field penetrating into the tachocline 
is consistent whenever r\ ^ 10 9 cm 2 /s. 
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1. Introduction 

Helioseismic inversions of the solar internal rotation during the past 
decade invariably showed that the surface-like latitudinal differential 
rotation, pervading the convective zone, changes to a near-rigid rota- 
tion in the radiative zone. The change takes place in a very thin layer 
known as the tachocline (e.g. Kosovichev, 1996). Known properties of 
the tachocline were recently reviewed by Corbard et al. (2001). While 
for a long time only upper limits were available for the thickness of the 
tachocline, a more exact determination of the thickness was recently 
made by several groups (Corbard et al., 1998; Corbard et al., 1999; 
Charbonneau et al., 1999; Basu and Antia, 2001). Their findings can 
be summarized as follows. Below the equator the tachocline is cen- 
tered around at R = 0.691 ± 0.003 solar radii, and its full thickness is 
w = 0.04 ±0.014 solar radii. 1 Compared with the internal radius of the 
(adiabatically stratified) convective zone (r^ cz = 0.71336 ± 0.00002), 

By "full thickness" here we mean the radial interval over which the horizontal 
differential rotation is reduced by a factor of 100. Different authors use different 
definitions of w, which explains most of the variation among various published values. 
The e-folding height of differential rotation, i.e. the scale height of the tachocline is 
then H = w/ In 100 <J 0.01 R Q . 
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this implies that the tachocline lies directly beneath the convective 
zone. There seems to be significant evidence for a slightly prolate 
form of the tachocline (R = 0.71 ± 0.003 solar radii at a latitude of 
60°) and marginal evidence for a thicker tachocline at high latitudes 
(w = 0.05 ± 0.005 solar radii at 60°). This seems to indicate that at 
higher latitudes up to half of the tachocline lies in the adiabatically 
stratified convective zone. 

A thorough understanding of the physics of the tachocline is crucial 
for understanding the solar dynamo for several reasons (see Petrovay, 
2000) . Firstly, the shear due to differential rotation is generally thought 
to be responsible for the generation of strong toroidal fields from poloidal 
fields. This shear is undoubtedly far stronger in the tachocline than 
anywhere else in the Sun. Second, both linear and nonlinear stability 
analyses of toroidal flux tubes lying at the bottom of the convective 
zone show that the storage of these tubes for times comparable to the 
solar cycle is only possible below the unstably stratified part of the 
convective zone, in a layer that crudely coincides with the tachocline. 
Indeed, as flux emergence calculations now provide ample evidence that 
solar active regions originate from the buoyant instability of 10 5 G flux 
tubes lying in the stably stratified layers below the convective zone 
proper (Moreno- Insertis, 1994), it is hard to evade the conclusion that 
strong magnetic fields oscillating with the dynamo period of 22 years 
must be present in (at least part of) the tachocline. 

A magnetic field oscillating with a circular frequency u> C yc = 
P = 22 years is known to penetrate a conductive medium only down 
to a skin depth of 

H skhl = (27/Kye) 1 / 2 (1) 

(cf. Garaud, 1999). Using a molecular value for the magnetic diffusivity 
T), this turns out to be very small (order of a few kilometers) in the solar 
case, apparently suggesting that the oscillating dynamo field cannot 
penetrate very deep into the tachocline region. Note, however, that 
in order to store a magnetic flux of order 10 23 Mx in the form of a 
toroidal field of 10 5 G in an active belt of width ~ 10 5 km, the storage 
region must clearly have a thickness of at least several megameters. 
(Or possibly more, taking into account the strong magnetic flux loss 
from emerging loops — cf. Petrovay and Moreno-Insertis, 1997; Dorch, 
private communication.) This is only compatible with equation (1) for 
r) ^ 10 9 cm 2 /s. 

A turbulent magnetic diffusivity of this order of magnitude is not 
implausible, given our present lack of detailed information about the 
physical conditions in the tachocline. The value required is still several 
orders of magnitude below the diffusivity in the solar convective zone 
(~ 10 13 cm 2 /s). The turbulence responsible for this diffusivity may be 
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generated either by non-adiabatic overshooting convection or by MHD 
instabilities of the tachocline itself. Note that several recent analyses 
have treated the problem of the stability of the tachocline under various 
approximations (e.g. Gilman and Dikpati, 2000; see review by Gilman, 
2000). The results are not conclusive yet, but it is clearly quite possible 
that, once all relevant (three-dimensional, nonlinear, MHD) effects are 
taken into account, the tachocline will prove to be unstable and there- 
fore capable of maintaining a certain level of turbulence. Note that the 
remaining slight inconsistencies between the standard and seismic solar 
models also seem to indicate some extra (probably turbulent) mixing in 
the tachocline layer (Gough, 2000). On the other hand, this turbulent 
diffusivity should not extend below a depth of a few times 10 Mm to 
avoid an overdepletion of lithium. Such a shallow depth of the turbulent 
layer is consistent with the seismic constraints on tachocline thickness, 
quoted above. 

These considerations prompt us to consider the structure of a tur- 
bulent tachocline pervaded by an oscillatory magnetic field. For sim- 
plicity, the poloidal field will be treated as given, in the form of a 
simple oscillating field of characteristic strength B p , of dipolar latitude- 
dependence, penetrating below the convective zone to a shallow depth 
of about 30 Mm (0.04 Rq). We will find that the observed tachocline 
thickness is reproduced for suitable pairs of the parameter values (B p , 
rf). This effect would offer a straightforward explanation for the thinness 
of the tachoclinc. Indeed, this thinness implies a horizontal transfer 
of angular momentum that is much more effective than the vertical 
transport. In our model, the horizontal transport is due to Maxwell 
stresses in the strong oscillatory magnetic field. Alternative mechanisms 
proposed for the horizontal transport include a strongly horizontally 
anisotropic turbulence (Spiegel and Zahn, 1992), non-diffusive hydro- 
dynamical momentum fluxes (Forgacs-Dajka and Petrovay, 2001) as 
well as Maxwell stresses due to a weak permanent internal magnetic 
field in the radiative interior (MacGregor and Charbonneau, 1999; 
Riidiger and Kitchatinov, 1997; Garaud, 2001). The plausibility of the 
processes implied by the the hydro dynamical scenarios is, however, 
dubious. As pointed out by Canuto (1998), the extremely strong hor- 
izontal anisotropy necessary to limit the tachocline to within 4 % of 
the solar radius has never been observed in nature or in laboratory, 
while non-diffusive fluxes also require an unrealistically high amplitude 
to do the trick. Internal magnetic fields are certainly able to confine a 
non-turbulent tachocline (and they will probably be needed to explain 
the lack of radial differential rotation in the solar interior anyway). 
Nevertheless, on basis of the arguments outlined above we feel that the 
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alternative concept of a turbulent tachocline pervaded and confined by 
a dynamo-generated field is worth considering. 

2. Estimates 

Let us first regard the following model problem. Consider a plane par- 
allel layer of incompressible fluid of density p, where the viscosity v and 
the magnetic diffusivity r\ are taken to be constant. At z = where z is 
the vertical coordinate (corresponding to depth in the solar application 
we have in mind) a periodic horizontal shearing flow is imposed in the 
y direction: 

v y0 = v o cos(kx) (2) 

(so that x will correspond to heliographic latitude, while y to the lon- 
gitude). We assume a two-dimensional flow pattern [d y = 0) and v x = 
v z = (no "meridional flow" ) . An oscillatory horizontal "poloidal" field 
is prescribed in the x direction as 

B x = B p cos(wi) (3) 

The evolution of the azimuthal components of the velocity and the 
magnetic field is then described by the corresponding components of 
the equations of motion and induction, respectively. Introducing v = v y 
and using Alfven speed units for the magnetic field 

V p = B^Ttp)- 1 ' 2 b = B y (47T P y 1 / 2 , (4) 

these can be written as 

dtv = V p cos(ujt)d x b — uV 2 v (5) 

dtb = V p cos(Lot)d x v — r]V 2 b (6) 
Solutions may be sought in the form 

v = v(x, z) + v'(x, z)f(u)t) (7) 

b = b'(x,z)f{ut + (j)) (8) 

where / is a 2-7r-periodic function of zero mean and of amplitude 0{l). 
(a denotes time average of a, while a' = a — a.) 

It will be of interest to consider the (temporal) average of equation 

(5): 

= V p cos{ujt)f(ujt + 4>)d x b' - uV 2 v (9) 

Subtracting this from equation (5) yields 

<V = Vp[cos{ut)f{ut + </>)]' d x b - uV 2 v' (10) 
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Figure 1. Magnetic field strength in Alfven speed units V p necessary to confine the 
tachocline to thickness H as a function of turbulent magnetic diffusivity r\: (a) for 
vjr\ — 1 with different values of H; (b) for H = O.Ol-R© with different values of 
the Prandtl number. (Note that, by coincidence, in the solar tachocline B p ~ V p to 
order of magnitude, in CGS units.) 



For an estimate, we will suppose cos(ut) f (u>t + 0) £ C(i) (i.e. 
no "conspiracy" between the phases, a rather natural assumption). 
As H <C R we may approximate V 2 ~ H~ 2 . Estimating the other 
derivatives as dt ~ oj and d x ~ i? -1 , (9) yields 



(11) 



A similar order-of-magnitude estimate of the terms in equation (10) 
yields 

(w + zV#V ~ V^&'/i? (12) 
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while from (6) we find in a similar manner 

ubf ~ (V p + v')V p /R + r^'/tf 2 



(13) 



From the last three order-of-magnitude relations one can work out with 
some algebra 



Equation (14) then tells us the field amplitude (in Alfven units) V p 
needed to confine the tachocline to a thickness H with given values 
of the diffusivities and of the oscillation period. If the diffusivities are 
of turbulent origin one expects v/rj ~ 1 and equation (14) can be 
used to plot V p as a function of r/ for different values of H, with uj = 
27r/22yrs (Fig. la). Figure lb illustrates the sensitivity of this result 
to our assumption about the magnetic Prandtl number v/r]. 

Of course, the (Vp, rj) pairs that reproduce the observed relation 
H ~ 5 Mm are also subject to the condition H 5s -ff s kin> otherwise the 
assumption of an oscillatory field pervading the tachocline would not 
be consistent. Using equation (1), this implies that only the regime 
r] ^ 10 9 cm 2 /s should be considered. 

Figure 1 thus suggests that an oscillatory poloidal field of about a 
thousand gauss is able to confine the tachocline to its observed thick- 
ness for rj ~ 10 10 cm 2 /s; for higher diffusivities a somewhat stronger 
field is needed. This figure can then serve as our guide in choosing the 
parameters of the more realistic spherical models considered in the next 
section. 



3.1. Equations 

In this section we focus on the numerical solution in the case of spherical 
geometry, using the realistic solar stratification. The time evolution of 
the velocity field v and magnetic field B are governed by the Navier- 
Stokes and induction equations: 




uR 2 uj (1 + r]/ujH 2 )(l + v/ojH 2 ) 
~H* 1 + 2u/ojH 2 



(14) 



3. Numerical Solution 



1 / B 2 \ 
$v + (v V)v = -W--V [P+— + 




Airp 



1 



(B-V)B + -V-r, (15) 



P 



and 



3tB = Vx(vxB)-Vx(i)Vx B), 



(16) 
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where r is the viscous stress tensor, V is the gravitational potential 
and p is pressure. These equations are supplemented by the constraints 
of mass and magnetic flux conservation 

V • (pv) = 0, (17) 

V-B = 0. (18) 

In our model we use the azimuthal component of the equations as- 
suming axial symmetry and we ignore the meridional flow. Then the 
magnetic and velocity fields can be written as 



B = -^— a d e A, l —d r A, B , (19) 

v = rsmOu)(r,0,t)e <t) , (20) 

where the usual spherical coordinates are used, A is the poloidal field 
potential, B is the toroidal field, lo is the angular velocity and is the 
azimuthal unit vector (cf. Rudiger and Kitchatinov, 1997). In order to 
present more transparent equations we write the poloidal field potential 
in the following form: 

A = a(r, t) sin 2 9. (21) 

The components of the viscous stress tensor appearing in the azimuthal 
component of the Navier-Stokes equation are 

re<t> = r^e = pv^- dg (4^) , (22) 
r V sin 9 J 

T<f> r = T r(t> = pVT 8 r (^y^j , (23) 

where v is the viscosity. 

Thus, the equations, including the effects of diffusion, toroidal field 
production by the differential rotation and the Lorentz force, read 

d t uJ = (d r v + 4- + -d r p) d r LO + vd 2 r uj + 3 "° OS d e uj + ^d 2 e uo (24) 
V r p J r 2 sm6 r z 

^ a cos 6 d r a ^ ^ ^ / a cos 9 d r acos9 \^ 

2tt pr 3 sin 9 r Airpr 3 \27rpr 4 sin0 47rpr 3 sin#/ ' 

2a sin 9 cos 9 d r asm 2 9 
OtB = o r uj Oquj (25) 

+ (*L + drV ) drB + r,d 2 B + ^- a d B + \d 2 B 

| ( d r y r] cos 2 9 _ Jl \ B 
\ r r 2 sin 2 9 r 2 
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Under the assumption of no meridional circulation used here, these 
equations remain unchanged when written in a reference frame rotating 
with an angular frequency Q. 

3.2. Boundary conditions 

The computational domain for the present calculations consists of just 
the upper part of the radiative interior, between radii r- m and r\^ cz . 
We use the same boundary conditions for uj(r,6,t) as those in Elliott 
(1997). We suppose that the rotation rate at the base of convection 
zone can be described with the same expression as in the upper part of 
the convection zone. In accordance with the observations of the GONG 
network, the following expression is given for f2b cz : 

^1 = 456 - 72 cos 2 9 - 42 cos 4 9 nHz. (26) 

This is used to give the outer boundary condition on u, 

0, + uj = S7 bcz at r = r hcz . (27) 

is chosen as the rotation rate in the radiative interior below the 
tachocline. On the basis of helioseismic measurements, this value is 
equal to the rotation rate of the convection zone at a latitude of about 
30°, corresponding to Q/2ir « 437 nHz. 

The second boundary condition oj = is imposed at the inner edge 
of our domain r; n . 

Our boundary conditions for the toroidal field are simply 

B = at r = rbcz and r = r- m . (28) 

3.3. Initial conditions 

The initial conditions chosen for all calculations are 



'JCZ 



u(r, 0, t = 0) = f^bcz ~ ^ at r = 
u(r,6,t = 0) = at r < r hcz (29) 

B(r,6,t = 0) = 0. 



3.4. Numerical method 

We used a time relaxation method with a finite difference scheme first 
order accurate in time to solve the equations. A uniformly spaced grid, 
with spacings Ar and A# is set up with equal numbers of points in the 
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Figure 2. The poloidal field configuration. 



r and 6 directions, r is chosen to vary from r- m = 4.2 x 10 u cm to rb cz , 
and 9 varies from to ir/2. 

The stability of this explicit time evolution scheme is determined by 
the condition 



where 



At < Min {At diff , Ai magdiff , At M axwcii} (30) 

rj \V 2 Ar 



Ar 2 Ar 2 

Atdiff = — At mag diflf = — AtMaxwell = 



2v n 



2% 



Here, B p is the amplitude of the poloidal field, defined here as 



Bp = Max 



4a 2 (d r a) 2 



bcz 



bcz 



(31) 



(32) 



Our calculations are based on a more recent version of the solar 
model of Guenther et al. (1992). 



4. Numerical Results and Discussion 

In this section we examine the influence of an oscillatory magnetic 
field on the radial spreading of the differential rotation on the basis of 
our numerical results. In our calculations the poloidal magnetic field is 
assumed to be time independent in amplitude, oscillating in time and 
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a priori known. For the function a of the poloidal field potential the 
following expression is used: 

a = A)COs(u,t)% (L^JElY r > nm (33) 

2 V r bcz r J 

a = r im > r, 

where ^4o fixes the field amplitude and n m is the depth of the pene- 
tration of the poloidal magnetic field into the radiative interior. As the 
skin effect should limit the penetration of the oscillatory field below 
the turbulent tachocline, we set r- im = 30 Mm. Figure 2 illustrates the 
poloidal field geometry. Note that the radial coordinate r shown on the 
ordinates in our figures has its zero point reset to r^ cz , i.e. r in the 
figures corresponds to r — r-^ cz . Thus, the negative r values correspond 
to the layers below the convective zone. 

The diffusive timescale over which the solution should relax to a 
very nearly periodic behaviour is 

T = (r bcz - nn) 2 A?mm (34) 

where ?y m i n is the lowest value of rj in the domain, i.e. its value taken at 
r m . Physically we would expect this value to be close to the molecular 
viscosity. Using this value in the computations would, however, lead to 
a prohibitively high number of timesteps to relaxation. (Clearly, the 
runtime of our computations must be chosen to well exceed r to reach 
relaxation.) 

Therefore, we first consider the simpler case 7] = 10 10 cm 2 /s = const., 
corresponding to r ~ 120 years. The results for this case are shown in 
Figure 3-5 for different amplitudes of the poloidal magnetic field, after 
relaxation. In accordance with the results of Section 2 it is found that 
a kilogauss poloidal field (peak amplitude 2400 G) is able to confine 
the tachocline to its observed thickness. 

Figure 6 presents the distribution of the azimuthal component of 
the Lorentz acceleration 

h = ^(B • V) B(f> (35) 

in a meridional section of the tachocline, showing polar acceleration and 
equatorial deceleration. The time variation of the differential rotation 
and the toroidal magnetic field during a cycle after relaxation is shown 
in Fig. 7. Note the significant variation with cycle phase. 

Finally, in the calculation presented in Figure 8-9 the diffusivity is 
allowed to vary with radius as 

rj = rjo exp ( _^ b ° z ) • (36) 
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Figure 3. Spreading of the differential rotation into the radiative interior for 
B p — 1603 G. Left-hand panel: contours of the time-average of the angular rota- 
tion rate u>(r,d,t) under one dynamo period. Equidistant contour levels are shown, 
separated by intervals of 100 nHz, starting from towards both non-negative (solid) 
and negative (dashed) values. Right-hand panel: differential rotation amplitude Alo 
(defined as the difference of maximal and minimal ui values in a horizontal surface) 
as a function of radius. 




This variation, while quantitatively much milder than expected physi- 
cally, allows us to consider the effects of a variable diffusivity without 
the necessity of prohibitively long integration times. The influence of 
the downwards decreasing diffusivity on the solution is found to be 
small. 

Now we turn to the comparison of the detailed spatio-temporal 
behaviour of our solutions to other models and helioseismic measure- 
ments. As a caveat we should stress that, in contrast to the above 
conclusions about tachocline confinement, these details of our models 
are unlikely to represent the situation in the real Sun faithfully, as the 
assumption of a simple dipole geometry (i.e. the use of a "standing 
wave" instead of a travelling dynamo wave) is clearly not realistic. Yet 
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Figure 6. Contours of the time-averaged azimuthal component of the Lorentz 
acceleration fL(r,6,t) for the case in Fig. 5. Contour levels correspond to 
log|/i| = (—6, —5.5, —5, —4.5, —4), solid contours standing for /l > 0, dashed 
contours for ft, < 0. 

it may be interesting from a theoretical point of view to draw some 
parallels with other work. In particular, the polar "cone" of rotational 
deceleration discernible in Figs. 3-7 is reminiscent of the similar feature 
seen in models with a steady internal field (cf. Fig. 2C of MacGregor 
and Charbonneau, 1999). This feature is, however, likely to disappear 
or change completely in a model with a poloidal field prescribed as a 
travelling dynamo wave. 

In Fig. 10 we present the timedatitude diagram of angular velocity 
fluctuations during one cycle. What is remarkable in this plot is not 
so much the weak poleward migration (which is sensitive to model 
details, and for other forms of A(r, 9, t) completely different migration 
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Elapsed time inside of a cycle: 11.72 year 
Poloidal field strength: -2350.50 Gauss 




Figure 7. Snapshot of the solution at four cycle phases for B p — 2405 G. t = cor- 
responds to poloidal field maximum. Left-hand panels: contours of angular rotation 
rate, as in Fig. 3. Right-hand panels: contours of toroidal magnetic field strength. 
Equidistant contour levels are shown, separated by intervals of 2000 G. 
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Figure 8. Radial profile of the diffusivity used for Figure 9. 
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Figure 9. The same as in Figure 5, but the diffusivity is varied in radius. 



patterns may result) but the fact that no other periodicity than the 
imposed 11/22 years is visible. In particular, there is no trace of the 
controversial 1.3 year periodicity recently claimed by several authors 
(Corbard et al, 2001 and references therein). 

Computer animations of the time evolution illustrated in Fig. 7 and 
other cases may be found on the accompanying CD-ROM, or they can 
be downloaded from the following internet address: 
http : / / astro . elte . hu/kutat/ sol/mgconf /mgconf e . html. 



5. Conclusion 

We have demonstrated that the tachocline can be confined to its ob- 
served thickness by a dipole-like poloidal magnetic field of about one 
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Figure 10. Time-latitude diagram of fluctuations of the angular velocity uj around 
its temporal average at a given latitude for the case in Fig. 5, at two different levels 
in the tachocline. Top: r - r bcz = -3.87 Mm; Bottom: r - r bcz = -27.09 Mm 

kilogauss, penetrating below the convective zone and oscillating with a 
period of 22 years, if the tachocline region is turbulent with a diffusivity 
of rj ~ 10 10 cm 2 /s. Our estimates presented in Section 2 suggest that a 
similar confinement may be produced for other pairs of the parameter 
values (Bp, rj). For the numerical computations we assumed a magnetic 
Prandtl number vjr]= 1; Fig. 1(b) gives an impression of the influence 
of using another v/rj value on the quantitative results 

A crucial assumption in these calculations is that the oscillatory 
poloidal magnetic field does penetrate below the convective zone to 
depths comparable to the tachocline thickness. This remains an as- 
sumption here, the penetration being prescribed a priori by setting 
the parameter r; m . While this should be confirmed later in more gen- 
eral calculations including the evolution equation for the poloidal field, 
formula (1) does suggest that the assumption is consistent whenever 
rj ^ 10 9 cm 2 /s. 

This explains the necessity of the high turbulent diffusivity in the 
tachocline, used in the present models. The high diffusivity, in turn, 
is the principal reason why the field strength necessary for tachocline 
confinement is so much higher here than in the case of models with a 
steady internal remnant magnetic field (MacGregor and Charbonneau, 
1999; Riidiger and Kitchatinov, 1997; Garaud, 2001): as the turbulent 
magnetic Prandtl number is not expected to differ greatly (i.e. by more 
than an order of magnitude) from unity, viscous angular momentum 
transfer is also quite effective. 
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Another effect that may contribute to the high field strength needed 
is that here the field is fully anchored in the differentially rotating over- 
lying envelope: since B r is nonzero at the upper boundary, the magnetic 
field applies stresses vertically as well as horizontally. In ideal MHD, for 
a stationary state, magnetic stresses are known to impose u> = const, 
along field lines (Ferraro's theorem). For finite but low diffusivities, 
field line anchoring can still make a great difference in the resulting 
rotational profile, as demonstrated by MacGregor and Charbonneau 
(1999). The fact that periodic oscillations of uj around its temporal 
mean in our model (Fig. 10) are strongest near the poles, where the field 
is nearly radial, shows that, despite the high diffusivity, this effect is to 
some extent also significant in our models. Thus, field line anchoring 
may be responsible for the fact that the field strengths required in the 
numerical models exceed somewhat the analytic estimates of Section 2, 
based on the assumption of a horizontal field. 

The ~ 1000 G poloidal field strength assumed here is nonetheless not 
implausible. Indeed, transport equilibrium models of the poloidal field 
inside the convective zone (Petrovay and Szakaly, 1999) show that the 
characteristic strength of the poloidal field in the bulk of the convective 
zone is order of 10 G; flux conservation then sets a lower limit of a few 
hundred gauss for the poloidal field strength in the region below the 
convective zone where the field lines close. This is quite compatible 
with peak values of 1-2000 G. 

Our neglect of meridional circulation may also seem worrying to 
some, given that Spiegel and Zahn (1992) found that the Eddington- 
Sweet circulation transports angular momentum much more efficiently 
than viscosity, and so it is the primary mechanism responsible for the 
inward spreading of the tachocline. This is, however, not expected to be 
the case for a turbulent tachocline. Owing to the strongly subadiabatic 
stratification below the convective zone, the timescale of any merid- 
ional circulation cannot be shorter than the relevant thermal diffusive 
timescale (to allow moving fluid elements to get rid of their buoyancy) . 
In the present case the relevant heat conductivity is the turbulent one, 
so for turbulent Prandtl numbers not very different from unity, the 
shortest possible timescale for meridional circulation is just compara- 
ble to the viscous timescale. Thus, while meridional circulation may 
possibly indeed have a significant effect on tachocline structure, in our 
view it is not likely to dominate over turbulent viscosity. 

It is clear that the models presented in this paper are still far from 
being a realistic representation of the solar tachocline. In particular 

- In order to limit the number of free parameters, we made no at- 
tempt here to simulate the latitudinal migration of the magnetic 
field during the cycle. 
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- The diffusivity, if generated by instabilities of the tachocline, should 
in some way be related to (e.g.) the local shear, instead of being 
arbitrarily specified. 

- The poloidal field evolution should also be consistently calculated 
in a more complete model. 

- Meridional circulation may also be taken into account. 

And the list could be continued. Despite these simplifications, we be- 
lieve that our model convincingly demonstrates the feasibility of tachocline 
confinement by a dynamo field. Possible generalizations of these mod- 
els can be made according to the remarks listed above. Work in this 
direction is in progress. 
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MAGNETIC FIELD 



E. Forgacs-Dajka and K. Petrovay 

Eotvos University, Dept. of Astronomy, Budapest, Pf. 32, H-1518 Hungary 
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Abstract. Helioseismic measurements indicate that the solar tachocline is very 
\^ ' thin, its full thickness not exceeding 4% of the solar radius. The mechanism that 

, inhibits differential rotation to propagate from the convective zone to deeper into 

the radiative zone is not known, though several propositions have been made. In 
£N) ' this paper we demonstrate by numerical models and analytic estimates that the 

tachocline can be confined to its observed thickness by a poloidal magnetic field 
Bp of about one kilogauss, penetrating below the convective zone and oscillating 
with a period of 22 years, if the tachocline region is turbulent with a diffusivity 
of 77 ~ 10 10 cm 2 /s (for a turbulent magnetic Prandtl number of unity). We also 
show that a similar confinement may be produced for other pairs of the parameter 
, values (Bp, rf). The assumption of the dynamo field penetrating into the tachocline 

^ • is consistent whenever r\ <; 10 9 cm 2 /s. 

CO . 

i Keywords: Sun: interior, MHD, tachocline 



CN 



O 

o 
23 



o 

C3 



X 



1. Introduction 



O Helioseismic inversions of the solar internal rotation during the past 

decade invariably showed that the surface-like latitudinal differential 
^ ' rotation, pervading the convective zone, changes to a near-rigid rota- 

tion in the radiative zone. The change takes place in a very thin layer 
known as the tachocline (e.g. Kosovichev, 1996). Known properties of 
the tachocline were recently reviewed by Corbard et al. (2001). While 
for a long time only upper limits were available for the thickness of the 



tachocline, a more exact determination of the thickness was recently 
made by several groups (Corbard et al., 1998; Corbard et al., 1999; 
Charbonneau et al, 1999; Basu and Antia, 2001). Their findings can 
be summarized as follows. Below the equator the tachocline is cen- 
tered around at R = 0.691 ± 0.003 solar radii, and its full thickness is 
w = 0.04 ±0.014 solar radii. 1 Compared with the internal radius of the 
(adiabatically stratified) convective zone (r^ cz = 0.71336 ± 0.00002), 



By "full thickness" here we mean the radial interval over which the horizontal 
differential rotation is reduced by a factor of 100. Different authors use different 
definitions of w, which explains most of the variation among various published values. 
The e-folding height of differential rotation, i.e. the scale height of the tachocline is 
then H = w/ In 100 <, 0.01 R Q . 

© 2008 Kluwer Academic Publishers. Printed in the Netherlands. 
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this implies that the tachocline lies directly beneath the convective 
zone. There seems to be significant evidence for a slightly prolate 
form of the tachocline (R = 0.71 ± 0.003 solar radii at a latitude of 
60°) and marginal evidence for a thicker tachocline at high latitudes 
(w = 0.05 ± 0.005 solar radii at 60°). This seems to indicate that at 
higher latitudes up to half of the tachocline lies in the adiabatically 
stratified convective zone. 

A thorough understanding of the physics of the tachocline is crucial 
for understanding the solar dynamo for several reasons (see Petrovay, 
2000) . Firstly, the shear due to differential rotation is generally thought 
to be responsible for the generation of strong toroidal fields from poloidal 
fields. This shear is undoubtedly far stronger in the tachocline than 
anywhere else in the Sun. Second, both linear and nonlinear stability 
analyses of toroidal flux tubes lying at the bottom of the convective 
zone show that the storage of these tubes for times comparable to the 
solar cycle is only possible below the unstably stratified part of the 
convective zone, in a layer that crudely coincides with the tachocline. 
Indeed, as flux emergence calculations now provide ample evidence that 
solar active regions originate from the buoyant instability of 10 5 G flux 
tubes lying in the stably stratified layers below the convective zone 
proper (Moreno- Insertis, 1994), it is hard to evade the conclusion that 
strong magnetic fields oscillating with the dynamo period of 22 years 
must be present in (at least part of) the tachocline. 

A magnetic field oscillating with a circular frequency u> C yc = 
P = 22 years is known to penetrate a conductive medium only down 
to a skin depth of 

H skhl = (27/Kye) 1 / 2 (1) 

(cf. Garaud, 1999). Using a molecular value for the magnetic diffusivity 
T), this turns out to be very small (order of a few kilometers) in the solar 
case, apparently suggesting that the oscillating dynamo field cannot 
penetrate very deep into the tachocline region. Note, however, that 
in order to store a magnetic flux of order 10 23 Mx in the form of a 
toroidal field of 10 5 G in an active belt of width ~ 10 5 km, the storage 
region must clearly have a thickness of at least several megameters. 
(Or possibly more, taking into account the strong magnetic flux loss 
from emerging loops — cf. Petrovay and Moreno-Insertis, 1997; Dorch, 
private communication.) This is only compatible with equation (1) for 
r) ^ 10 9 cm 2 /s. 

A turbulent magnetic diffusivity of this order of magnitude is not 
implausible, given our present lack of detailed information about the 
physical conditions in the tachocline. The value required is still several 
orders of magnitude below the diffusivity in the solar convective zone 
(~ 10 13 cm 2 /s). The turbulence responsible for this diffusivity may be 
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generated either by non-adiabatic overshooting convection or by MHD 
instabilities of the tachocline itself. Note that several recent analyses 
have treated the problem of the stability of the tachocline under various 
approximations (e.g. Gilman and Dikpati, 2000; see review by Gilman, 
2000). The results are not conclusive yet, but it is clearly quite possible 
that, once all relevant (three-dimensional, nonlinear, MHD) effects are 
taken into account, the tachocline will prove to be unstable and there- 
fore capable of maintaining a certain level of turbulence. Note that the 
remaining slight inconsistencies between the standard and seismic solar 
models also seem to indicate some extra (probably turbulent) mixing in 
the tachocline layer (Gough, 2000). On the other hand, this turbulent 
diffusivity should not extend below a depth of a few times 10 Mm to 
avoid an overdepletion of lithium. Such a shallow depth of the turbulent 
layer is consistent with the seismic constraints on tachocline thickness, 
quoted above. 

These considerations prompt us to consider the structure of a tur- 
bulent tachocline pervaded by an oscillatory magnetic field. For sim- 
plicity, the poloidal field will be treated as given, in the form of a 
simple oscillating field of characteristic strength B p , of dipolar latitude- 
dependence, penetrating below the convective zone to a shallow depth 
of about 30 Mm (0.04 Rq). We will find that the observed tachocline 
thickness is reproduced for suitable pairs of the parameter values (B p , 
rf). This effect would offer a straightforward explanation for the thinness 
of the tachoclinc. Indeed, this thinness implies a horizontal transfer 
of angular momentum that is much more effective than the vertical 
transport. In our model, the horizontal transport is due to Maxwell 
stresses in the strong oscillatory magnetic field. Alternative mechanisms 
proposed for the horizontal transport include a strongly horizontally 
anisotropic turbulence (Spiegel and Zahn, 1992), non-diffusive hydro- 
dynamical momentum fluxes (Forgacs-Dajka and Petrovay, 2001) as 
well as Maxwell stresses due to a weak permanent internal magnetic 
field in the radiative interior (MacGregor and Charbonneau, 1999; 
Riidiger and Kitchatinov, 1997; Garaud, 2001). The plausibility of the 
processes implied by the the hydro dynamical scenarios is, however, 
dubious. As pointed out by Canuto (1998), the extremely strong hor- 
izontal anisotropy necessary to limit the tachocline to within 4 % of 
the solar radius has never been observed in nature or in laboratory, 
while non-diffusive fluxes also require an unrealistically high amplitude 
to do the trick. Internal magnetic fields are certainly able to confine a 
non-turbulent tachocline (and they will probably be needed to explain 
the lack of radial differential rotation in the solar interior anyway). 
Nevertheless, on basis of the arguments outlined above we feel that the 
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alternative concept of a turbulent tachocline pervaded and confined by 
a dynamo-generated field is worth considering. 

2. Estimates 

Let us first regard the following model problem. Consider a plane par- 
allel layer of incompressible fluid of density p, where the viscosity v and 
the magnetic diffusivity r\ are taken to be constant. At z = where z is 
the vertical coordinate (corresponding to depth in the solar application 
we have in mind) a periodic horizontal shearing flow is imposed in the 
y direction: 

v y0 = v o cos(kx) (2) 

(so that x will correspond to heliographic latitude, while y to the lon- 
gitude). We assume a two-dimensional flow pattern [d y = 0) and v x = 
v z = (no "meridional flow" ) . An oscillatory horizontal "poloidal" field 
is prescribed in the x direction as 

B x = B p cos(wi) (3) 

The evolution of the azimuthal components of the velocity and the 
magnetic field is then described by the corresponding components of 
the equations of motion and induction, respectively. Introducing v = v y 
and using Alfven speed units for the magnetic field 

V p = B^Ttp)- 1 ' 2 b = B y (47T P y 1 / 2 , (4) 

these can be written as 

dtv = V p cos(ujt)d x b — uV 2 v (5) 

dtb = V p cos(Lot)d x v — r]V 2 b (6) 
Solutions may be sought in the form 

v = v(x, z) + v'(x, z)f(u)t) (7) 

b = b'(x,z)f{ut + (j)) (8) 

where / is a 2-7r-periodic function of zero mean and of amplitude 0{l). 
(a denotes time average of a, while a' = a — a.) 

It will be of interest to consider the (temporal) average of equation 

(5): 

= V p cos{ujt)f(ujt + 4>)d x b' - uV 2 v (9) 

Subtracting this from equation (5) yields 

<V = Vp[cos{ut)f{ut + </>)]' d x b - uV 2 v' (10) 
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Figure 1. Magnetic field strength in Alfven speed units V p necessary to confine the 
tachocline to thickness H as a function of turbulent magnetic diffusivity r\: (a) for 
vjr\ — 1 with different values of H; (b) for H = O.Ol-R© with different values of 
the Prandtl number. (Note that, by coincidence, in the solar tachocline B p ~ V p to 
order of magnitude, in CGS units.) 



For an estimate, we will suppose cos(ut) f (u>t + 0) £ C(i) (i.e. 
no "conspiracy" between the phases, a rather natural assumption). 
As H <C R we may approximate V 2 ~ H~ 2 . Estimating the other 
derivatives as dt ~ oj and d x ~ i? -1 , (9) yields 



(11) 



A similar order-of-magnitude estimate of the terms in equation (10) 
yields 

(w + zV#V ~ V^&'/i? (12) 
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while from (6) we find in a similar manner 

ubf ~ (V p + v')V p /R + r^'/tf 2 



(13) 



From the last three order-of-magnitude relations one can work out with 
some algebra 



Equation (14) then tells us the field amplitude (in Alfven units) V p 
needed to confine the tachocline to a thickness H with given values 
of the diffusivities and of the oscillation period. If the diffusivities are 
of turbulent origin one expects v/rj ~ 1 and equation (14) can be 
used to plot V p as a function of r/ for different values of H, with uj = 
27r/22yrs (Fig. la). Figure lb illustrates the sensitivity of this result 
to our assumption about the magnetic Prandtl number v/r]. 

Of course, the (Vp, rj) pairs that reproduce the observed relation 
H ~ 5 Mm are also subject to the condition H 5s -ff s kin> otherwise the 
assumption of an oscillatory field pervading the tachocline would not 
be consistent. Using equation (1), this implies that only the regime 
r] ^ 10 9 cm 2 /s should be considered. 

Figure 1 thus suggests that an oscillatory poloidal field of about a 
thousand gauss is able to confine the tachocline to its observed thick- 
ness for rj ~ 10 10 cm 2 /s; for higher diffusivities a somewhat stronger 
field is needed. This figure can then serve as our guide in choosing the 
parameters of the more realistic spherical models considered in the next 
section. 



3.1. Equations 

In this section we focus on the numerical solution in the case of spherical 
geometry, using the realistic solar stratification. The time evolution of 
the velocity field v and magnetic field B are governed by the Navier- 
Stokes and induction equations: 




uR 2 uj (1 + r]/ujH 2 )(l + v/ojH 2 ) 
~H* 1 + 2u/ojH 2 



(14) 



3. Numerical Solution 



1 / B 2 \ 
$v + (v V)v = -W--V [P+— + 




Airp 



1 



(B-V)B + -V-r, (15) 



P 



and 



3tB = Vx(vxB)-Vx(i)Vx B), 



(16) 
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where r is the viscous stress tensor, V is the gravitational potential 
and p is pressure. These equations are supplemented by the constraints 
of mass and magnetic flux conservation 

V • (pv) = 0, (17) 

V-B = 0. (18) 

In our model we use the azimuthal component of the equations as- 
suming axial symmetry and we ignore the meridional flow. Then the 
magnetic and velocity fields can be written as 



B = -^— a d e A, l —d r A, B , (19) 

v = rsmOu)(r,0,t)e <t) , (20) 

where the usual spherical coordinates are used, A is the poloidal field 
potential, B is the toroidal field, lo is the angular velocity and is the 
azimuthal unit vector (cf. Rudiger and Kitchatinov, 1997). In order to 
present more transparent equations we write the poloidal field potential 
in the following form: 

A = a(r, t) sin 2 9. (21) 

The components of the viscous stress tensor appearing in the azimuthal 
component of the Navier-Stokes equation are 

re<t> = r^e = pv^- dg (4^) , (22) 
r V sin 9 J 

T<f> r = T r(t> = pVT 8 r (^y^j , (23) 

where v is the viscosity. 

Thus, the equations, including the effects of diffusion, toroidal field 
production by the differential rotation and the Lorentz force, read 

d t uJ = (d r v + 4- + -d r p) d r LO + vd 2 r uj + 3 "° OS d e uj + ^d 2 e uo (24) 
V r p J r 2 sm6 r z 

^ a cos 6 d r a ^ ^ ^ / a cos 9 d r acos9 \^ 

2tt pr 3 sin 9 r Airpr 3 \27rpr 4 sin0 47rpr 3 sin#/ ' 

2a sin 9 cos 9 d r asm 2 9 
OtB = o r uj Oquj (25) 

+ (*L + drV ) drB + r,d 2 B + ^- a d B + \d 2 B 

| ( d r y r] cos 2 9 _ Jl \ B 
\ r r 2 sin 2 9 r 2 
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Under the assumption of no meridional circulation used here, these 
equations remain unchanged when written in a reference frame rotating 
with an angular frequency Q. 

3.2. Boundary conditions 

The computational domain for the present calculations consists of just 
the upper part of the radiative interior, between radii r- m and r\^ cz . 
We use the same boundary conditions for uj(r,6,t) as those in Elliott 
(1997). We suppose that the rotation rate at the base of convection 
zone can be described with the same expression as in the upper part of 
the convection zone. In accordance with the observations of the GONG 
network, the following expression is given for f2b cz : 

^1 = 456 - 72 cos 2 9 - 42 cos 4 9 nHz. (26) 

This is used to give the outer boundary condition on u, 

0, + uj = S7 bcz at r = r hcz . (27) 

is chosen as the rotation rate in the radiative interior below the 
tachocline. On the basis of helioseismic measurements, this value is 
equal to the rotation rate of the convection zone at a latitude of about 
30°, corresponding to Q/2ir « 437 nHz. 

The second boundary condition oj = is imposed at the inner edge 
of our domain r; n . 

Our boundary conditions for the toroidal field are simply 

B = at r = rbcz and r = r- m . (28) 

3.3. Initial conditions 

The initial conditions chosen for all calculations are 



'JCZ 



u(r, 0, t = 0) = f^bcz ~ ^ at r = 
u(r,6,t = 0) = at r < r hcz (29) 

B(r,6,t = 0) = 0. 



3.4. Numerical method 

We used a time relaxation method with a finite difference scheme first 
order accurate in time to solve the equations. A uniformly spaced grid, 
with spacings Ar and A# is set up with equal numbers of points in the 



mgconf.tex; 1/02/2008; 18:25; p. 8 



Tachoclinc Confinement by an Oscillatory Magnetic Field 



9 




-20 



-30 : 



-40 ; 



-50 : 



-60 E 



r im = -30 Mm 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 
[rod] 



Figure 2. The poloidal field configuration. 



r and 6 directions, r is chosen to vary from r- m = 4.2 x 10 u cm to rb cz , 
and 9 varies from to ir/2. 

The stability of this explicit time evolution scheme is determined by 
the condition 



where 



At < Min {At diff , Ai magdiff , At M axwcii} (30) 

rj \V 2 Ar 



Ar 2 Ar 2 

Atdiff = — At mag diflf = — AtMaxwell = 



2v n 



2% 



Here, B p is the amplitude of the poloidal field, defined here as 



Bp = Max 



4a 2 (d r a) 2 



bcz 



bcz 



(31) 



(32) 



Our calculations are based on a more recent version of the solar 
model of Guenther et al. (1992). 



4. Numerical Results and Discussion 

In this section we examine the influence of an oscillatory magnetic 
field on the radial spreading of the differential rotation on the basis of 
our numerical results. In our calculations the poloidal magnetic field is 
assumed to be time independent in amplitude, oscillating in time and 
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a priori known. For the function a of the poloidal field potential the 
following expression is used: 

a = A)COs(u,t)% (L^JElY r > nm (33) 

2 V r bcz r J 

a = r im > r, 

where ^4o fixes the field amplitude and n m is the depth of the pene- 
tration of the poloidal magnetic field into the radiative interior. As the 
skin effect should limit the penetration of the oscillatory field below 
the turbulent tachocline, we set r- im = 30 Mm. Figure 2 illustrates the 
poloidal field geometry. Note that the radial coordinate r shown on the 
ordinates in our figures has its zero point reset to r^ cz , i.e. r in the 
figures corresponds to r — r-^ cz . Thus, the negative r values correspond 
to the layers below the convective zone. 

The diffusive timescale over which the solution should relax to a 
very nearly periodic behaviour is 

T = (r bcz - nn) 2 A?mm (34) 

where ?y m i n is the lowest value of rj in the domain, i.e. its value taken at 
r m . Physically we would expect this value to be close to the molecular 
viscosity. Using this value in the computations would, however, lead to 
a prohibitively high number of timesteps to relaxation. (Clearly, the 
runtime of our computations must be chosen to well exceed r to reach 
relaxation.) 

Therefore, we first consider the simpler case 7] = 10 10 cm 2 /s = const., 
corresponding to r ~ 120 years. The results for this case are shown in 
Figure 3-5 for different amplitudes of the poloidal magnetic field, after 
relaxation. In accordance with the results of Section 2 it is found that 
a kilogauss poloidal field (peak amplitude 2400 G) is able to confine 
the tachocline to its observed thickness. 

Figure 6 presents the distribution of the azimuthal component of 
the Lorentz acceleration 

h = ^(B • V) B(f> (35) 

in a meridional section of the tachocline, showing polar acceleration and 
equatorial deceleration. The time variation of the differential rotation 
and the toroidal magnetic field during a cycle after relaxation is shown 
in Fig. 7. Note the significant variation with cycle phase. 

Finally, in the calculation presented in Figure 8-9 the diffusivity is 
allowed to vary with radius as 

rj = rjo exp ( _^ b ° z ) • (36) 
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Figure 3. Spreading of the differential rotation into the radiative interior for 
B p — 1603 G. Left-hand panel: contours of the time-average of the angular rota- 
tion rate u>(r,d,t) under one dynamo period. Equidistant contour levels are shown, 
separated by intervals of 100 nHz, starting from towards both non-negative (solid) 
and negative (dashed) values. Right-hand panel: differential rotation amplitude Alo 
(defined as the difference of maximal and minimal ui values in a horizontal surface) 
as a function of radius. 




This variation, while quantitatively much milder than expected physi- 
cally, allows us to consider the effects of a variable diffusivity without 
the necessity of prohibitively long integration times. The influence of 
the downwards decreasing diffusivity on the solution is found to be 
small. 

Now we turn to the comparison of the detailed spatio-temporal 
behaviour of our solutions to other models and helioseismic measure- 
ments. As a caveat we should stress that, in contrast to the above 
conclusions about tachocline confinement, these details of our models 
are unlikely to represent the situation in the real Sun faithfully, as the 
assumption of a simple dipole geometry (i.e. the use of a "standing 
wave" instead of a travelling dynamo wave) is clearly not realistic. Yet 
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Figure 6. Contours of the time-averaged azimuthal component of the Lorentz 
acceleration fL(r,6,t) for the case in Fig. 5. Contour levels correspond to 
log|/i| = (—6, —5.5, —5, —4.5, —4), solid contours standing for /l > 0, dashed 
contours for ft, < 0. 

it may be interesting from a theoretical point of view to draw some 
parallels with other work. In particular, the polar "cone" of rotational 
deceleration discernible in Figs. 3-7 is reminiscent of the similar feature 
seen in models with a steady internal field (cf. Fig. 2C of MacGregor 
and Charbonneau, 1999). This feature is, however, likely to disappear 
or change completely in a model with a poloidal field prescribed as a 
travelling dynamo wave. 

In Fig. 10 we present the timedatitude diagram of angular velocity 
fluctuations during one cycle. What is remarkable in this plot is not 
so much the weak poleward migration (which is sensitive to model 
details, and for other forms of A(r, 9, t) completely different migration 
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Elapsed time inside of a cycle: 11.72 year 
Poloidal field strength: -2350.50 Gauss 




Figure 7. Snapshot of the solution at four cycle phases for B p — 2405 G. t = cor- 
responds to poloidal field maximum. Left-hand panels: contours of angular rotation 
rate, as in Fig. 3. Right-hand panels: contours of toroidal magnetic field strength. 
Equidistant contour levels are shown, separated by intervals of 2000 G. 
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Figure 8. Radial profile of the diffusivity used for Figure 9. 
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Figure 9. The same as in Figure 5, but the diffusivity is varied in radius. 



patterns may result) but the fact that no other periodicity than the 
imposed 11/22 years is visible. In particular, there is no trace of the 
controversial 1.3 year periodicity recently claimed by several authors 
(Corbard et al, 2001 and references therein). 

Computer animations of the time evolution illustrated in Fig. 7 and 
other cases may be found on the accompanying CD-ROM, or they can 
be downloaded from the following internet address: 
http : / / astro . elte . hu/kutat/ sol/mgconf /mgconf e . html. 



5. Conclusion 

We have demonstrated that the tachocline can be confined to its ob- 
served thickness by a dipole-like poloidal magnetic field of about one 
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Figure 10. Time-latitude diagram of fluctuations of the angular velocity uj around 
its temporal average at a given latitude for the case in Fig. 5, at two different levels 
in the tachocline. Top: r - r bcz = -3.87 Mm; Bottom: r - r bcz = -27.09 Mm 

kilogauss, penetrating below the convective zone and oscillating with a 
period of 22 years, if the tachocline region is turbulent with a diffusivity 
of rj ~ 10 10 cm 2 /s. Our estimates presented in Section 2 suggest that a 
similar confinement may be produced for other pairs of the parameter 
values (Bp, rj). For the numerical computations we assumed a magnetic 
Prandtl number vjr]= 1; Fig. 1(b) gives an impression of the influence 
of using another v/rj value on the quantitative results 

A crucial assumption in these calculations is that the oscillatory 
poloidal magnetic field does penetrate below the convective zone to 
depths comparable to the tachocline thickness. This remains an as- 
sumption here, the penetration being prescribed a priori by setting 
the parameter r; m . While this should be confirmed later in more gen- 
eral calculations including the evolution equation for the poloidal field, 
formula (1) does suggest that the assumption is consistent whenever 
rj ^ 10 9 cm 2 /s. 

This explains the necessity of the high turbulent diffusivity in the 
tachocline, used in the present models. The high diffusivity, in turn, 
is the principal reason why the field strength necessary for tachocline 
confinement is so much higher here than in the case of models with a 
steady internal remnant magnetic field (MacGregor and Charbonneau, 
1999; Riidiger and Kitchatinov, 1997; Garaud, 2001): as the turbulent 
magnetic Prandtl number is not expected to differ greatly (i.e. by more 
than an order of magnitude) from unity, viscous angular momentum 
transfer is also quite effective. 
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Another effect that may contribute to the high field strength needed 
is that here the field is fully anchored in the differentially rotating over- 
lying envelope: since B r is nonzero at the upper boundary, the magnetic 
field applies stresses vertically as well as horizontally. In ideal MHD, for 
a stationary state, magnetic stresses are known to impose u> = const, 
along field lines (Ferraro's theorem). For finite but low diffusivities, 
field line anchoring can still make a great difference in the resulting 
rotational profile, as demonstrated by MacGregor and Charbonneau 
(1999). The fact that periodic oscillations of uj around its temporal 
mean in our model (Fig. 10) are strongest near the poles, where the field 
is nearly radial, shows that, despite the high diffusivity, this effect is to 
some extent also significant in our models. Thus, field line anchoring 
may be responsible for the fact that the field strengths required in the 
numerical models exceed somewhat the analytic estimates of Section 2, 
based on the assumption of a horizontal field. 

The ~ 1000 G poloidal field strength assumed here is nonetheless not 
implausible. Indeed, transport equilibrium models of the poloidal field 
inside the convective zone (Petrovay and Szakaly, 1999) show that the 
characteristic strength of the poloidal field in the bulk of the convective 
zone is order of 10 G; flux conservation then sets a lower limit of a few 
hundred gauss for the poloidal field strength in the region below the 
convective zone where the field lines close. This is quite compatible 
with peak values of 1-2000 G. 

Our neglect of meridional circulation may also seem worrying to 
some, given that Spiegel and Zahn (1992) found that the Eddington- 
Sweet circulation transports angular momentum much more efficiently 
than viscosity, and so it is the primary mechanism responsible for the 
inward spreading of the tachocline. This is, however, not expected to be 
the case for a turbulent tachocline. Owing to the strongly subadiabatic 
stratification below the convective zone, the timescale of any merid- 
ional circulation cannot be shorter than the relevant thermal diffusive 
timescale (to allow moving fluid elements to get rid of their buoyancy) . 
In the present case the relevant heat conductivity is the turbulent one, 
so for turbulent Prandtl numbers not very different from unity, the 
shortest possible timescale for meridional circulation is just compara- 
ble to the viscous timescale. Thus, while meridional circulation may 
possibly indeed have a significant effect on tachocline structure, in our 
view it is not likely to dominate over turbulent viscosity. 

It is clear that the models presented in this paper are still far from 
being a realistic representation of the solar tachocline. In particular 

- In order to limit the number of free parameters, we made no at- 
tempt here to simulate the latitudinal migration of the magnetic 
field during the cycle. 
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- The diffusivity, if generated by instabilities of the tachocline, should 
in some way be related to (e.g.) the local shear, instead of being 
arbitrarily specified. 

- The poloidal field evolution should also be consistently calculated 
in a more complete model. 

- Meridional circulation may also be taken into account. 

And the list could be continued. Despite these simplifications, we be- 
lieve that our model convincingly demonstrates the feasibility of tachocline 
confinement by a dynamo field. Possible generalizations of these mod- 
els can be made according to the remarks listed above. Work in this 
direction is in progress. 
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